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[i] A theory of cross-coupled flow equations in unsaturated soils is necessary to predict 
(1) electroosmotic flow with application to electroremediation and agriculture, (2) the 
electroseismic and the seismoelectric effects to develop new geophysical methods to 
characterize the vadose zone, and (3) the streaming current, which can be used to 
investigate remotely ground water flow in unsaturated conditions in the capillary water 
regime. To develop such a theory, the cross-coupled generalized Darcy and Ohm 
constitutive equations of transport are extended to unsaturated conditions. This model 
accounts for inertial effects and for the polarization of porous materials. Rather than using 
the zeta potential, like in conventional theories for the saturated case, the key parameter 
used here is the quasi-static volumetric charge density of the pore space, which can be 
directly computed from the quasi-static permeability. The apparent permeability entering 
Darcy' s law is also frequency dependent with a critical relaxation time that is, in turn, 
dependent on saturation. A decrease of saturation increases the associated relaxation 
frequency. The final form of the equations couples the Maxwell equations and a simplified 
form of two-fluid phases Biot theory accounting for water saturation. A generalized 
expression of the Richard equation is derived, accounting for the effect of the vibration of 
the skeleton during the passage of seismic waves and the electrical field. A new expression 
is obtained for the effective stress tensor. The model is tested against experimental data 
regarding the saturation and frequency dependence of the streaming potential coupling 
coefficient. The model is also adapted for two-phase flow conditions and a numerical 
application is shown for water flooding of a nonaqueous phase liquid (NAPL, oil) 
contaminated aquifer. Seismoelectric conversions are mostly taking place at the NAPL 
(oil)/water encroachment front and can be therefore used to remotely track the position of 
this front. This is not the case for other geophysical methods. 

Citation: Revil, A., and H. Mahardika (2013), Coupled hydromechanical and electromagnetic disturbances in unsaturated porous 
materials, Water Resour. Res., 49, doi: 10.1002/wrcr.20092. 



1. Introduction 

[2] Revil et al. [20 1 1 ] recently developed a new set of 
constitutive equations to model cross-coupled transport 
phenomena in porous media. This model was however re- 
stricted to fully water saturated materials and for quasi- 
static flow problems (inertial effects were neglected). The 
motivation of the present paper is to extend the macro- 
scopic generalized (cross-coupled) Darcy and Ohm laws to 
unsaturated porous materials and to account for dynamic 
effects (harmonic pressure or harmonic electrical fields and 
inertial terms). The final equations couple Biot's theory for 
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unsaturated porous media to the Maxwell equations. The 
coupling is electrokinetic in nature. 

[3] There are broad applications of such a model for the 
assessment of water resources and remediation. For 
instance, the record of the electrical field of electrokinetic 
nature can be used to image nonintrusively ground water 
flow and to get access easily to the parameters characteriz- 
ing the capillary pressure curve and the relative permeabil- 
ity [Jougnot et al., 2012; Mboh et al., 2012]. With a cross- 
coupled flow theory, it is also possible to model electroos- 
mosis, which can be used to move solutes and nonaqueous 
phase liquids (NAPLs)/dense NAPLs (DNAPLs) in the 
vadose zone and aquifers for remediation purposes \Acar 
and Alshawabkeh, 1993; Bruell et al., 1992; Han et al., 
2004]. This method can also be used to dewater clayey 
soils in civil engineering [Bjerrum, 1967]. Finally, electro- 
osmosis can be used also to move moisture up to the roots 
of plants [Huweg et al., 2010], the required electrical field 
can be produced through the use of photovoltaics [Kamel, 
1994]. 

[4] The theory developed herein can be used to predict 
the electroseismic (electric-to- seismic conversion) effect in 
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a broad range of frequencies (from 0.1 Hz to 100 kHz). The 
electroseismic coupling is related to the generation of seis- 
mic waves (hydromechanical disturbances) in response to a 
nonstationary electrical field [Reppert and Morgan, 2002; 
Thompson, 2005; Thompson et ah, 2005]. While this effect 
has been mainly used so far for oil-related problems, it may 
be used to map NAPLs and DNAPLs as well [Hinz, 201 1]. 
Alternatively, the seismoelectric (seismic-to-electric, SE) 
effect corresponds to the generation of electromagnetic dis- 
turbances from seismic waves [e.g., Ivanov, 1939; Frenkel, 
1944; Dupuis and Butler, 2006; Dupuis et al., 2009; 
Jardani et al., 2010]. Hydromechanical disturbances (e.g., 
Haines jumps in unsaturated flow conditions, see Haas and 
Revil [2009]) can also be responsible for the generation of 
electromagnetic disturbances. These disturbances can in 
turn be remotely monitored by electromagnetic methods. 
To our knowledge, this is the first theory able to predict the 
electroseismic and SE effects in unsaturated conditions 
with applications to vadose zone hydrogeology. Dupuis 



et al. [2007] has demonstrated the feasibility to measure SE 
conversions for vadose zone applications. 

2. Generalized Constitutive Equations 
2.1. Assumptions 

[5] We consider below an isotropic porous body with 
connected pores. The surface of the minerals in contact 
with water is negatively charged at near-neutral pH values. 
There is therefore an excess of electrical charges in the 
pore space. This excess of charges occurs in the electrical 
double layer (Figure 1). This electrical double layer is 
actually made up of two layers: (1) a layer of (counter) ions 
sorbed directly onto the mineral surface and (2) a diffuse 
layer. In the following, we will use the subscript a to describe 
the properties of air, and subscripts w and s for the water and 
solid phases, respectively. Water is assumed to be the wetting 
phase. The term skeleton (or frame) is used to describe the 
assemblage of grains alone without fluids in the pores. 
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Figure 1. Double layer associated with clays, (a) TOT structure of clay minerals, where T stands for 
the tetrahedral layer (mainly Si) and O stands for the octohedral layer (mainly Al). (b) The surface of 
clay minerals in contact with water is charged because of both the existence of amphoteric sites on the 
edges of the clay crystal and basal negative sites associated with isomorphic substitutions in the crystal- 
line framework of the clay minerals in the T layer and O layer, (c) The mineral charge is compensated 
by counterions (M + ) and coions (A~) forming a double layer. This double layer comprises a layer of 
sorbed counterions (the Stern layer) and a diffuse layer where the Coulombic interactions between the 
charged mineral surface and the coions and counterions prevail. Note that a similar double layer exists at 
the surface of silica due to the reactivity of silanol sites with water. The parameter/ denotes the fraction 
of the countercharge located in the Stern layer (partition coefficient). 
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[6] Another set of assumptions used below pertain to the 
capillary pressure curve. Hysteretic behavior will be 
neglected, and therefore the porous material will thus be 
characterized by a unique set of hydraulic functions. In sec- 
tions 1-5, we also work with the Richards model that 
assumes that the air pressure is constant and equal to the 
atmospheric pressure. This implies in turn that the air phase 
is infinitely mobile and connected to the atmosphere (in 
section 7, we will consider the case of two-phases flow 
using the approach of Rubino and Holliger [2012]). The 
capillary pressure p c (in Pascal) is defined as the pressure 
of the nonwetting phase minus the pressure of the wetting 
phase [e.g., Bear and Verruijt, 1987], 



Table 1. Nomenclature of the Nonmechanical Material Properties 



-Pa ~P* 



(1) 



where p a and p w denote the average air and water pressures 
(in Pascal), respectively. The capillary head (in meter) is 
defined as 



Pc^ 



- Pa 



(2) 



In unsaturated flow conditions, the gradient of the capillary 
head is given by 



W = — Vp w . 



(3) 



This assumption is used to avoid dealing with the flow of 
the air phase. In unsaturated conditions, the capillary pres- 
sure is positive, the capillary head (suction) is negative, 
and the pressure of the water phase is smaller than the 
atmospheric pressure. The total head includes the gravity 
force and is defined by ip + z, where z denotes the elevation. 
In the following, the model will be restricted to the capillary 
regime (saturation comprised between the irreducible water 
saturation and full saturation). Extending the present model 
below irreducible saturation would require a unified water 
capillary curve model (including a water sorption isotherm) 
and a description of film flow along the surface of the grains. 
This is outside the scope of the present paper. 

[7] Finally, attenuation of the seismic waves associated 
with squirt-flow dissipation mechanisms will be neglected 
despite the fact that this mechanism is known to control the 
attenuation of seismic waves in the frequency band usually 
used for seismic field investigations (see Rubino and Hol- 
liger [2012] for a recent pore-scale modeling of this effect 
and Dvorkin et al. [1994] and Carcione [2007] for a 
description of the relevant physics). 

2.2. Generalized Darcy's Law 

[8] In saturated conditions, the (averaged) filtration dis- 
placement is defined as [e.g., Morency and Tromp, 2008], 



(4) 



where u and u M . denote the averaged displacement of the solid 
and water phases, respectively, and <fi denotes the connected 
porosity. A nomenclature of the material properties and me- 
chanical properties used below is provided in Tables 1 and 2, 
respectively. The disturbances are considered to be harmonic 
of the form exp(— iujt), where i = s/—l denotes the pure 
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imaginary number, u> = 2nf denotes the angular frequency (in 
rad s _1 ), and / the frequency (in Hertz). In saturated condi- 
tions, the Darcy velocity is defined as the time derivative of 
the filtration displacement when the skeleton is moving. In 
water- saturated conditions, the generalized Darcy's law is 
given by [e.g., Jardani et al., 2010], 



- — (Vp + p r ii + pfFyv 

Vf V 



(5) 



Table 2. Nomenclature of the Mechanical Material Properties 
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where 0 = d@/dt denotes the time derivative of the pa- 
rameter 0 and t represents time, k 0 denotes the low-fre- 
quency permeability of the porous material (in m 2 ), F 
denotes the electrical formation factor, which is the ratio of 
the bulk tortuosity of the pore space to the connected poros- 
ity (see for instance, Pride [1994]), Ff denotes the body 
force applied to the pore water phase (in N m -3 , e.g., the 
gravitational body force or the electrical force acting on the 
excess of electrical charges of the pore water). The forma- 
tion factor is related to the porosity by the first Archie's 
law F = <\T"' ', where m is called the cementation exponent 
[Archie, 1942]. To keep the notations as light as possible, 
we will not distinguish below if the variables are expressed 
in the time domain or in the frequency domain. It is easy to 
recognize if the equations are written in the frequency do- 
main or in the time domain, the switch from one domain to 
the other being done by a Fourier transform or its associ- 
ated inverse Fourier transform. 

[9] In unsaturated conditions, the filtration displacement 
and the mass density of the water phase are given by, 



[11] The complex permeability is used to describe the 
effect of the inertial force in Darcy's law [e.g., Pride, 
1994]. In the time domain, a nonlinear Darcy's law can be 
similarly defined in which the apparent permeability 
depends on the Reynolds number (see Boleve et al. [2007] 
for the development of the equations and an experimental 
check). Alternatively, this equation could be written to be 
consistent with the Forchheimer equation. In the For- 
chheimer equation, the flux is decomposed in several terms 
while in the equations derived in Boleve et al. [2007], this 
is the right-hand side of the Darcy equation that is devel- 
oped as a nonlinear function of the fluid pressure (nonlinear 
Darcy equation). Aulisa et al. [2009] derived a straightfor- 
ward approach to go from one formulation to the other. The 
relaxation time characterizes the transition between the 
viscous laminar flow regime and the inertial laminar flow 
regime in the Navier-Stokes equation formulated in the fre- 
quency domain. The critical frequency associated with this 
relaxation time is given by, 



w w 



(6) 



A 



2nT k 2nk r k 0 p w F 



(12) 



Pf 



:t 



v)Pg + S w p w , 



(7) 



respectively, s w denotes the water saturation (s w = 1 at satu- 
ration). The mass density of the gas phase can be neglected, 
and therefore p f ~ s w p w . From equation (4), the porosity 
can be replaced by s w <f> when dealing with unsaturated con- 
ditions (of course, terms in (1 — 0), dealing with the solid 
phase, remain unchanged). The Darcy velocity associated 
with the water phase is given by, 



Vw 



VPw + PwV" + P vt Sw(Fs w '')vf w - F„ 



(8) 



where k r denotes the relative permeability (dimensionless), 
r] w denotes the dynamic viscosity of the pore water (in Pa 
s), and p w is the pressure of the water phase (it will be 
replaced later by the suction head, see equation (3) above). 
The term Fs~" denotes the ratio of the bulk tortuosity of 
the water phase divided by the connected porosity, n 
denotes Archie's second or saturation Archie's exponent 
(n > 1, dimensionless, Archie [1942]). 

[10] From now, the constitutive equations are described 
in the frequency domain and w M , in the right-hand side of 
equation (8) can be replaced by — ium w and combined with 
the w w term found in the left-hand side of equation (8). 
Equation (8) can then be rewritten as, 

-ium w = - (Vp w - co 2 p w s n ,u - F w ) , (9) 
Vw 

where k* (u) is a complex- valued apparent permeability 
given by 



k r ko 



1 — ILOTk 

with the relaxation time t> given by, 

k r kop w F !_„ 

T k = s w ■ 

Vw 



(10) 



(11) 



Note that because n> 1 [Archie, 1942], n - 1 > 0. The 
measurement of this relaxation frequency can be used to 
estimate the (low-frequency) permeability of the material. 

[12] The two-flow regimes in porous media are sketched 
in Figure 2. The low-frequency flow regime (wt{ <C 1, 
/ "C /i) corresponds to the viscous laminar flow regime for 
which the flow in a cylindrical pore obeys Poiseuille's law 
(Figure 2c). High frequencies (or* 3> 1, / 3> fk) corre- 
sponds to the inertial laminar flow regime for which the 
pore water flow is described as being a potential-flow prob- 
lem. We consider here Reynolds number that are much 
below the values corresponding to the turbulent regime. 

[13] In poroelasticity, it is customary to define the fol- 
lowing two variables [e.g., Morency and Tromp, 2008; Jar- 
dani et al., 2010, and references therein], 



Vw 

k 0 



Pw = Pw F , 



(13) 



(14) 



where p w denotes an apparent pore water mass density. The 
relationship between the permeability and the water satura- 
tion can be expressed with the Brooks and Corey [1964] 
relationship, 



2+3A 

k — t A 



(15) 



From equation (11) and equations (13)— (15), we obtain, 



b 



(16) 



where g = 2/A + 4 — n. In Revil [2013, Appendix A], it is 
shown that (2 + 3 A) /A scales approximately as (n+2), so 
»«3. 
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a. Thick double layer b. Thin double layer 




c. Flow in the viscous laminar regime d. Flow in the inertial laminar regime 




Figure 2. Sketch of the charge distribution and flow regime in a pore. There are four end members to 
consider, depending on the pore size with respect to the double-layer thickness and on the frequency, (a) 
Thick double layer. The counterions of the diffuse layer are uniformly distributed in the pore space), (b) 
Thin double layer (the thickness of the diffuse layer is much smaller than the size of the pores), (c) Vis- 
cous laminar flow regime occurring at low frequencies, (d) Inertial laminar flow regime occurring at 
high frequencies. 



[14] In order to write a hydrodynamic equation coupled 
with the electrical field, the body force F w entering equa- 
tion (9) should be expressed by Coulomb's law, 

F„. = Q* v (uj)E, (17) 

where Qy(u) denotes the frequency-dependent excess of 
charge that can be dragged by the flow of the pore water 
through the pore space of the material (dynamic excess 
charge density of the pore space, see Appendix A) and E 
denotes the electrical field (in V m _1 ). The charge density 

Q v {oS) is frequency dependent because there is more 
charges dragged in the inertial laminar flow regime than in 
the viscous laminar flow regime in agreement with the 

model of Pride [1994]. In the following, the parameters Q v 

and Q v are the (dynamic or effective) volumetric charge 
density dragged in the low- frequency (wr t « 1) and high- 
frequency (ujTk» 1) regimes, respectively. Because the 
transition between low-frequency and high-frequency 
regimes is governed by the relaxation time t>, the following 



functional can be used to compute the effective charge den- 
sity as a function of the frequency (Appendix A) : 

Ea^) = W + {^~w) v^i^' (18) 

The low-frequency and high-frequency charge densities 

~ 0 ~ 00 

Q v and Q v can be found as follows. 
[15] (1) At low frequencies only a small fraction of the 

counterions of the diffuse layer are dragged by the flow of 

~ 0 ~ 00 

the pore water and therefore Q v « Q v . An expression to 

compute Q v directly from the low- frequency permeability 
&o is discussed further in section 5a below. 

[16] (2) At high frequencies, all the charge density exist- 
ing in the pores is uniformly dragged along the pore water 
flow, and therefore the charge density Q v is also equal to 
the volumetric charge density of the diffuse layer. An 
expression to compute Q v from the cation exchange 
capacity (CEC) is discussed further in section 5b below. 
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[17] Depending on the size of the electrical double layer 
with respect to the size of the pores, two cases need to be 
considered : 

[is] (1) In the thick double-layer approximation (see 

Figure 2a), Q v ~ Q v (all the counterions of the diffuse 
layer are dragged by the flow whatever the frequency), and 
therefore 



(19) 



contribution to the current density (source current density) 
is given by (see Appendix A), 



3 S = -iuQ r (u)wZ. 



(23) 



The mechanical contribution to the filtration displacement 
is given by the generalized Darcy's law derived in section 
2.2 above, 



(Vp» 



(24) 



[19] (2) In the thin double flayer approximation (see 
Figure 2b), one can expect Q v ^> Q v (see examples in 
Jougnot et al. [2012] and discussion in section 5 below) 
and therefore (see Appendix A) 



Q* r (ui,s w ) m QyS^y/l - iu>Tk{s w ). 



(20) 



[20] Introducing Coulomb's law, equation (17), into the 
Darcy equation, equation (9), yields, 



-iujv/ w = [Vp w - uj p w s w u) 



E. (21) 



Equation (21) shows the influence of three forcing terms on 
the Darcy velocity : (i) the displacement of the solid frame- 
work (through viscous drag at the pore water/solid inter- 
face), (ii) the pore fluid pressure gradient (in which a 
gravitational component can be added if needed), and (iii) 
the electrical field (which corresponds to a body force per 
unit charge density) through electroosmosis. Electroosmo- 
sis refers to the flow of the pore water in response to an 
electrical field. Physically, the positive excess of charge is 
responsible for a net flow in the direction of the electrical 
field through viscous drag of the pore water by the excess 
of electrical charge in the pore water. 

2.3. Generalized Ohm's Law 

[21] We investigate now the macroscopic electrical cur- 
rent density J. The first contribution is the conduction cur- 
rent density J e given by Ohm's law, 



J e = tr*ME, 



(22) 



where the conductivity a* (uj) denotes the complex conduc- 
tivity. For clayey materials, this conductivity has been 
modeled recently by Revil [2012, 2013]. The frequency de- 
pendence of the apparent conductivity is a direct conse- 
quence of the fact that the force applied to the charge 
carriers is controlled by an electrochemical potential, not 
just by the electrical field. Therefore, electromigration and 
diffusion are always coupled phenomena in porous media. 
The second contribution to the total current density corre- 
sponds to the advective drag of the excess of charge of the 
pore space by the flow of the pore water (contribution of 
advective nature). If the Darcy velocity associated with the 
poromechanical contribution is written as w'", the second 



[22] Equations (22)-(24) yield the following generalized 
Ohm's law, 



J = a* ME - 



k*(u)Q V (i0) 



(Vp w - ui 2 p w s w u) 



(25) 



[23] A model for the complex conductivity is now dis- 
cussed. The complex conductivity is written as, 



+ ia" = \a\exp(iip), 



(26) 



1 /2 

where | cr* | = (o / 2 + a" 2) ' denotes the magnitude of the 
conductivity and tp = atan(V' /V) denotes the phase lag 
between the electrical current and the resulting electrical 
field. For clayey materials, the complex conductivity is 
observed to be relatively independent on the frequency in 
the range 0.1 Hz-10 kHz [see Vinegar and Waxman, 1984; 
Revil, 2012, 2013]. The in-phase electrical conductivity is 
given as a function of the pore water conductivity a w by, 



's- 



(27) 



where F denotes the formation factor introduced above and 
a s denotes the surface conductivity. The surface conductiv- 
ity is given by the model developed by Revil [2012, 2013] : 



OS 



''w M 



" OC 

(+)Qv - 



A(6,m) 



where A(cf>, m) (~ 1 ) is defined by, 



A(4>, m) 



m(F- 1) 



3 VI -cj) 



(28) 



(29) 



(30) 



An alternative expression is given by Revil [2013]. 

[24] Equation (29) implies that the surface conductivity 
is controlled by the diffuse layer containing the fraction of 
counterions (1 —J) (Figure 1) characterized by an ionic mo- 
bility fi, + s equal to the mobility of the cations in the bulk 
pore water [Revil, 2012, 2013]. 

[25] Following Revil [2012, 2013], the quadrature con- 
ductivity can be expressed as, 



(3D 
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a" = 



'{+) 



f 
1-./ 



e 



(32) 



Cp(w, 



(39) 



where denotes the mobility of the counterions in the 
Stern layer (see value in Revil [2012, 2013]). These equa- 
tions provide a simple and accurate model to describe the 
complex conductivity of shaly sands and soils and general- 
ize the Vinegar and Waxman [1984] model. As noted by 
Vinegar and Waxman [1984] and Revil [2012], the fre- 
quency dependence of the quadrature conductivity is not 
explicit in equation (32). That said, for quasi-static condi- 
tions, the quadrature conductivity should go to zero under 
DC conditions, 



limer"(w) 



(33) 



The typical frequency below which the quadrature conduc- 
tivity becomes frequency dependent is typically smaller 
than 0.1 Hz [see Revil, 2012] and is therefore not relevant 
to the SE and SE problems. Note thus that the present 
approach includes a complex conductivity that agree with 
experimental data, while the frequency dependence of the 
complex conductivity in Pride [1994] is due to an electro- 
osmotic contribution that is expected to be negligible (see 
discussion in Marshall and Madden [1959]) and therefore 
cannot explain the observed quadrature conductivity of po- 
rous clayey materials (note that Pride [1994] stated very 
clearly in his paper that induced polarization was neglected 
in his approach). 

2.4. Coupled Constitutive Transport Equations 

[26] The two constitutive equations for the generalized 
Ohm and Darcy laws are written into the following matrix 
form: 



L*(lo,s w ) 

Vw 



where the cross-coupling term L* (w) is defined as, 

k*(u,s w )Q v(u,s w ) 



L*(u,s w ) 



(34) 



(35) 



The generalized streaming potential coupling coefficient is 
defined by the following equations in the quasi-static limit 
of Maxwell equations : 



lim Cp(u,s v ) 



(40) 



Equation (40) (developed by Revil and coworkers, see 
Revil et al. [2007] and Linde et al. [2007] for quasi-static 
streaming potentials) has been successfully tested by Mboh 
et al. [2012] through laboratory measurements. A check of 
this model will be provided below in section 6. 

[27] Similarly, the generalized electroosmotic coupling 
coefficient is defined in the quasi-static limit of the Max- 
well equations and when the skeleton is at rest and in ab- 
sence of pore fluid flow by the ratio between the gradient of 
pore fluid pressure divided by the gradient in electrical 
potential. For the thin-double layer case, this yields, 



I*(w,*„ 



>'/„■ 



C* 



lLOTk. 



(41) 



(42) 



(43) 



VxE = 0 = E = -Vtp, 



(36) 



[28] For the case corresponding to the thick double layer, 
the electroosmotic coupling coefficient is given by 

* —Q°y/ S w- Therefore, the electroosmotic coupling 
coefficient is simply a measure of the excess of charges 
that can be moved by the flow of the pore water and that 
this excess of charge increases with the frequency between 
the viscous laminar flow and the inertial laminar flow re- 
gime. It also increases with the decrease of the water satu- 
ration. Note however that in equation (34), this is L* (not 

Q v ) that controls the magnitude of the Darcy velocity, and 
therefore the intensity of the electroosmotic contribution to 
the Darcy velocity decreases with the saturation, which is 
qualitatively in agreement with the observations made by 
Aggour et al. [1994]. 

[29] In unsaturated flow conditions, it is customary to use 
the capillary head gradient V-0 = Vp w /p w g instead of the 
pore water pressure gradient (see section 2.1 above). Below 
the relaxation frequency separating the low-frequency vis- 
cous laminar flow regime from the high-frequency inertial 
laminar flow regime (true in the frequency range used for 
field applications), equation (34) can be rewritten in the time 
domain using the pressure head and including the gravita- 
tional field in the hydraulic driving force. This yields, 



difi\ L*(uj,. 
3?Vj=0;u=o o*^,, 

rj w a*(u!,s w ) 



(37) 



(38) 



More explicitly, in the thin-double layer approximation, 
the dynamic streaming potential coupling coefficient is 
given by, 



L* (s w ) 
k r (s w )k 0 



(44) 



where the quasi-static coupling coefficient L* is given by, 



k r (s w )k 0 Q v 



(45) 
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Equations (44) and (45) therefore provide a simple model 
for modeling the occurrence of streaming potential and 
electroosmosis in porous media. 

3. Description of the Hydromechanical Model 

3.1. Generalized Diffusion Equation for the Pore 
Water 

[30] The starting point is the following Biot constitutive 
equation in saturated conditions : 



-p = Cv V ■ u + M s V ■ w. 
Equation (46) is also often written as 

-p = C s e kk - M S C, 



(46) 



(47) 



where e kk = V ■ u = 5V/V (where V denotes the volume 
of the representative elementary volume) denotes the volu- 
metric strain of the porous body and ( = — V • w = 
4>(V • u — V • u„ ) denotes the linearized increment of fluid 
content [e.g., Lo et al., [2002]. The parameter £ represents 
the fractional volume of water flowing in or out of the rep- 
resentative volume of skeleton in response to an applied 
stress. The bulk moduli M s and C„ in saturated conditions, 
are defined as [e.g., Pride, 1994], 



(48) 



(1 -s w )/K a + s w /K w [Wood, 1955; Teja and Rice, 1981], 
where K a and K w represent the bulk moduli of air and 
water, respectively (K a = 0.145 MPa and K w = 2.25 GPa, 
see Lo et al. [2005])). 

[32] In unsaturated conditions, from equations (50) and 
(51), C= a w M and the Biot coefficient a w is given by 
where a denotes the Biot coefficient in saturated 
conditions. That said, there should be no exchange of water 
below the irreducible water saturation and therefore the 
correct scaling should be a w = s e a rather than a w = s w a, 
where s e denotes the reduced water saturation s e = 
(s w — s,.)/(l — s r ) where s r denotes the irreducible water 
saturation. In other words ( 1 /M) — > 0 as s w — > s r . 

[33] In addition to the poroelastic change corresponding 
to equation (5 1 ), there is also a change related to the mois- 
ture capacity. Generalizing equation (46) to unsaturated 
conditions yields therefore, 



1 d6 \ C 

M + dVj {P "~ P " )= M S/ - U + V - W - (53) 



where d9/dp w denotes the specific moisture capacity which 
is determined from the derivative of the capillary pressure 
with respect to the water content (in unsaturated flow, the 
air pressure is kept constant and close to saturation 
]xm(dQfdp v/ ) s _+ l — > 0). From section 2a, the filtration dis- 
placement of the water phase is given by 

-ium w = - (y Pw _ u 2 Pw s n ,u ~ F w ) . (54) 



1 



1+A 

K f 



(49) 



where a = 1 — Kf r /Ks denotes the Biot coefficient in the 
saturated state and A = (K f /(j>K s 2 ) [(1 - (j))K s ~ K fl ] . The 
Biot modulus M corresponds to the inverse of the poroelas- 
tic component of the specific storage and is defined as the 
increase of fluid (per unit volume of rock) as a result of an 
increase of pore pressure under constant volumetric strain. 
The following relationship applies C s = aM s . 

[31] To extend these equations to the unsaturated case, 
we apply the classical change of variables used in section 2 
above (p f s w p w , p ^ p w , 4>^ s w <f>, and w w w = 
s w (j)(u w — u)). This yields 



1 _ 1+A 

C ~ K f + K S A ' 



1 

M 



1+A 

K f 



(50) 



(51) 



where 9 = s n ,<f) denotes the water content (dimensionless) 
and where, in unsaturated conditions, the fluid increment is 
defined by [Berryman et al., 1988], 



-V • w w = <f>s w {V ■ u - V 



(52) 



The term A = (Kf/^Ks 1 ) [(1 - <f>)K s - K fr ] depends on 
the saturation because the compressibility of the fluid is 
given for instance by the Wood formula (\/Kf) = 



Therefore, the filtration displacement is given by, 

k*(io) 



(\7p w - uj 2 p w s w u - F w ) 



(55) 



Equations (53) and (55) yield, 



M <9p v 



(p a ~Pw) = a»V • u + 



V ■ 



(56) 



where the relationship C/M = a w has been used. Equation 
(56) is a nonlinear diffusion equation for the fluid pressure. 
For this to be obvious, the terms of this equation need to be 
reworked. Multiplying all the terms by (iujr] w /k*(uj)), sepa- 
rating the pressure terms in the left-hand side from the 
source term in the right-hand side and taking into consider- 
ation that the air pressure is constant (unsaturated flow 
assumption), the following nonlinear hydraulic diffusion 
equation is obtained, 



V ■ 



V ■ 



Vp w 



1 de\, 



k*(u) , 2 , 



zwa M ,V-u + V- ( — ^-F, 



(57) 
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Equation (57) can be written in the time domain using an 
inverse Fourier transform. Assuming that the permeability 
is given by its low-frequency asymptotic limit (which is 
correct below 10 kHz) and using Coulomb law plus the 
gravity force as body force (the frequency-dependent volu- 
metric charge density is also taken in its low frequency 
limit too) yield, 



V ■ 



-VP* 



i 



oo 



M dp„ 



k r k 0 p w s w 



a w V • u + V • 



k r k Q I g° 



E + P n .g 



(58) 



The origin of the three forcing terms on the right-hand side 
of this equation is now clearly established: the first term is 
related to the acceleration of the seismic wave acting on the 
skeleton of the material, the second term is due to the veloc- 
ity of the seismic wave, and the third term (at constant grav- 
ity acceleration) corresponds to the pore water flow 
associated with the electroosmotic forcing associated with 
the drag of the pore water by the electromigration of the 
excess of charge contained in the pore space of the material. 

[34] Another possibility is to write a generalized Rich- 
ards equation [Richards, 1931] showing the influence of 
the forcing terms in this equation (we assume again that the 
air pressure is constant). Starting with equation (58) and 
replacing the water pressure by the capillary head defined 
by ip — (p w — Pa)/p w g (in meter), the following Richards 
equation is obtained, 



K h — u 

g 



- a w V • u - V 



P v g «w 



(59) 



poroelastic term should be kept to have a formulation that 
remains consistent with the saturated state. The hydraulic 
conductivity is related to the relative permeability k r and 
Ko, the hydraulic conductivity at saturation. With the 
Brooks and Corey [1964] model, the porous material is sat- 
urated when the fluid pressure reaches the atmospheric 
pressure (-0 = 0 at the water table). The effective saturation, 
the specific moisture capacity, the relative permeability, 
and the water content are defined by 



96» 

b\p' 



(a b ip) , ip<l/a b , 

1, V > 

-Aa A (^-^)( ai V) (A+1) , V> 

4' > I/"*, 




{a b ip) 



■ + s e (4>-e r ), 



-(2+3A) 



ip < l/a b , 
ip > l/a b , 



ip < I/a*, 
ip > !/«*, 



(63) 



(64) 



(65) 



(66) 



respectively, where a b denotes the inverse of the capillary 
entry pressure related to the matric suction at which pore 
fluid begins to leave a drying porous material, A is called 
the pore size distribution index (a textural parameter), 9 r 
represents the residual water content (9 r = s r (p). Sometimes 
the residual water saturation is not accounted for and the 
capillary pressure curve and the relative permeability are 
written as, 



(67) 




Q^ + V-[-^(VV+l)] = -a w V-u + V 



A',, 



*wP w g 



Kh = kM^ = krKsj 



(60) 



(61) 



where a w = s e a = s e (l — Kfr/Ks) and where K = Kf r 
denotes the bulk modulus of the skeleton (drained bulk 
modulus), K h denotes the hydraulic conductivity (in m 
s _1 ), K s denotes the hydraulic conductivity at saturation, 
and C e denotes the specific storage term. This storage term 
is the sum of the specific moisture capacity (in m~ ) (also 
called the water capacity function) and the specific storage 
corresponding to the poroelastic deformation of the mate- 
rial. This yields, 



dip 



Pwg 

M 



(62) 



Usually in unsaturated conditions, the poroelastic term is 
much smaller than the specific moisture capacity, but the 



K 



Pc > Pe 
Pc < Pe 



(68) 



Because a w = s e a = s e (l — Kfr/Ks), when the water satu- 
ration reaches the irreducible water saturation, the two 
source terms on the right-hand side of equation (60) are 
null. Therefore, there is no possible excitation below the ir- 
reducible water saturation. In reality, this is not necessarily 
true, and the model could be completed by including film 
flow below the irreducible water saturation. 



3.2. Newton's Law and the Definition of the Effective 
Stress 

[35] The hydromechanical equations are defined in terms 
of an effective stress tensor. As explained in detail by Jar- 
dani et al. [2010] and Revil and Jardani [2010], there is a 
computational advantage in expressing the coupled hydro- 
mechanical problem in terms of the fluid pressure and dis- 
placement of the solid phase (four unknowns in total) 
rather than using the displacement of the solid and filtration 
displacement (six unknowns in total). 
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[36] In saturated conditions, Newton's law (which is a 
momentum conservation equation for the skeleton partially 
filled with its pore water) is written as, 



(69) 



V • T + F = -uf[pa + p f vi 



where T is the total stress tensor (positive normal stress for 
tension, see Detournay and Cheng [1993]) and F denotes 
the total body force applied to the porous material. In un- 
saturated conditions, Newton's law is written as, 



V • T + F : 



-uj 2 (pu + s w p w w w ), 



(70) 



where the filtration displacement of the pore water phase is 
given by, 



k*{uj) 



(Vp w - ui 2 p w s w u - F„ 



(71) 



[37] The last fundamental constitutive equation needed 
to complete the hydromechanical model in unsaturated 
conditions is a relationship between the total stress tensor 
(or the effective stress tensor) and the displacement of the 
solid phase and filtration displacement of the pore water 
phase. This equation is Hooke's law, which, in linear 
poroelasticity and for saturated conditions, is given by 



T = (A„V ■ u + CV ■ w)7 + G(Vu + Vu 7 



(80) 



where the infinitesimal deformation tensor is related to the 
displacement of the solid phase by e = (l/2)(Vu + Vu r ), 
G = Gfr denotes the shear modulus that is equal to the 
shear modulus of the skeleton (frame), and A„ = K u — 
(2/3)G denotes the Lame modulus in undrained conditions 
(K u denotes the undrained bulk modulus). In unsaturated 
conditions and accounting for the air pressure, equation 
(78) can be written as, 



Combining equations (70) and (71) yields, 



V • T + F : 



pa + s w p„ 



\7p w - u) p w s w a - F w ) 



V • T + F = -iJ-pla - lu 2 
where, 



k*{uS) , 



o , 2 k*UJ) 2 



(72) 
(73) 

(74) 



The effective stress in unsaturated conditions is taken as 
(this choice is justified below), 



T +pj + s e a(p w -p tt )I, 



(75) 



T+pJ = (A„V • u + CV • w w )7 + G(Vu + Vu 7 ) . (81) 

From equation (53) of section 3.1 above, the linearized in- 
crement of fluid content is given by, 



-V-w w =i(p w -^ a )+^V-u. 



Combining equations (81) and (82) yields, 



(82) 



T+Pa 



K u -~g\v ■u+c(-~( Pw -p a ) 



a„,V • u 



G(Vu + Vu 7 



T+pJ 



K u - a w C — -G ) V ■ u — a w (p w - p a ) 
G(Vu + Vu r ), 



(83) 



(84) 



which is both consistent with Bishop effective stress princi- 
ple in unsaturated conditions and the Biot stress principle 
in saturated conditions. The confining pressure and the 
effective confining pressure are defined as, 



P = Trace T, 

3 



-Trace r eff , 



(76) 



(77) 



respectively. This yields P e s = P — p a — s e a(p w —p a ) in 
agreement with the recent results by Lu et al. [2010]. Equa- 
tions (73) and (75) yield, 



_ = i k*{oj) 

V ■ T eS + F - u/s w p w ^^F„ 



iffftu + OaVp*, (78) 



where the hydromechanical coupling term 9 U is defined by 



2 k *{u) 
s w \a-uj p 



(79) 



where the following expression derived in section 3.1 
has been used for the Biot coefficient in unsaturated 
conditions, 



C 

— = a w = s e a. 
M 



(85) 



In addition, the bulk modulus is given by K = K u — a w C 
and the Lame Modulus is given by A = K — (2/3) G. Equa- 
tion (84) yields the following Hooke's law for the effective 
stress 

f + p a l + a w (p w - Pa)! = A(V ■ u)7 + G(Vu + Vu 7 ) , (86) 

T eS = A(V-u)7 + G(Vu + Vu r ), (87) 

where the effective stress is also given by equation (75). 
Note that the effective stress is only related to the deforma- 
tion of the skeleton of the porous material (by definition). 
This model generalizes therefore the effective stress con- 
cept developed recently by Lu et al. [2010]. 
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4. Description of Maxwell Equations 

[38] Pride [1994] volume averaged the local Maxwell 
equations to obtain a set of macroscopic Maxwell equations 
in the thin-double layer limit. The same equations were 
obtained by Revil and Linde [2006] for the thick double- 
layer case. The form of the four macroscopic Maxwell 
equations (Faraday's law of induction, Ampere's law, 
Gauss's law for magnetism, Gauss's law) is, 

VxE = -B, (88) 

VxH = J + D, (89) 

V ■ B = 0, (90) 



5. Determination of Charge Densities 

[40] The goal of this section is to provide a way to esti- 

- 0 - oc 

mate the two charge densities Q v and Q v used in the pre- 
vious sections. We first start by the low-frequency charge 
density and then we provide a model to estimate the high- 
frequency charge density. 

5.1. Determination of the Quasi-Static Excess Charge 
Density 

[41] For a fully water saturated material, the streaming 
potential coupling coefficient is defined as, 

C 0 = lim C(w, s w = 1) = . (96) 



V • D = p f , (91) 

where B = V x A denotes the magnetic field (in T), H is 
the auxiliary magnetic field (in A m~ ), and D is the current 
displacement vector (in C m -2 ), A is the magnetic (vector) 
potential, E = —V<p — A, where ip denotes the electric 
(scalar) potential (in V), and p f denotes the density of free 
charges These equations are completed by two electromag- 
netic constitutive equations : D = eE and B = fiH, where e 
is the permittivity of the porous body and u denotes its 
magnetic permeability. In the absence of magnetized 
grains, /i = /x 0 where denotes the magnetic permeability 
of free space. 

[39] When the harmonic electrical field is written as 
E = Eo exp(— iujt), and uj is the angular frequency with E 0 a 
constant electrical field magnitude and direction, the dis- 
placement current density vector is given by ij = D = 
— iuieE. The total current density 3 T entering Ampere's law, 



Equation (96) provides a way to estimate the charge density 

Q v from the measurements of the low-frequency streaming 
potential coupling coefficient, the low-frequency electrical 
conductivity, and permeability using, 

6;=%^°. (97) 

The estimate of the low-frequency volumetric charge 

density Q v is reported as a function of the permeability in 
Figure 3 for experiments performed with a broad range of 
porous materials at near-neutral pH values (pH 5-8). 

According to Jardani et al. [2007], Q v can be directly esti- 
mated from the quasi-static (saturated) permeability by (see 
Figure 3): 

logio Q°y = -9.23 - 0.82 log 10 k 0 . (98) 



is given by, 



V x H = Jt 



J + J,/, 



(92) 



(93) 



J r = a * (w)E + k MQL (y pw _ u? Pw s w \x) - iueK, (94) 

Vw 



Equation (98) therefore provides a way to estimate directly 
the volumetric charge density from the low-frequency per- 
meability, thereby reducing the number of material proper- 
ties to consider in the simulations. 

5.2. Determination of the High-Frequency Excess 
Charge Density 

[42] We now seek a way to estimate the high-frequency 
charge density Q v . At full saturation, the surface conduc- 
tivity is given by, 



J T = (a * (uj) - iLue)E + - Np w 
Vw 



■ ui p w s w u) 



(95) 



and an effective conductivity can be introduced such as 
°eff = a * ( w ) — iuJ£ ' where cr* ff = <r e ff + ;a>£ e ff is the 
effective or apparent complex conductivity and er e ff and 
e e ff are real scalars dependent upon frequency. These 
effective parameters are the parameters that are measured 
during an experiment in the laboratory or in the field and 
contain both electromigration and true dielectric polariza- 
tion effects. They are given by er e ff = a 1 and 
e e ff = o"/w — e. Another recently published paper is dedi- 
cated to the description of these parameters in the fre- 
quency range 1 mHz-1 GHz [Revil, 2013]. 



(99) 



Therefore, the determination of the surface conductivity 
and the formation factor (from the measurements of the 
conductivity of the porous material at different pore water 

conductivities) can be used to assess the value of . 

Therefore, the high-frequency excess charge density Q v 
can be estimated to the surface conductivity and the forma- 
tion factor by, 



a s F 



(100) 



Indeed, at high frequencies, all the charge density existing 
in the pores is uniformly dragged along the pore water 
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Figure 3. Quasi-static charge density Q° v (excess pore charge moveable by the quasi-static pore water 
flow) versus the quasi-static permeability for a broad collection of core samples and porous materials. 
This charge density is directly derived from the streaming potential coupling coefficient using equation 
(97). Data from Ahmad [1969], Boleve et al. [2007], Casagrande [1983], Friborg [1996], Jougnot et al. 
[2012], Jardani et al. [2007], Pengra et al. [1999], Revil et al. [2012, 2007], Sheffer [2007]; Revil et al. 
[2013], and Zhu and Toksoz [2012]. 



flow, and therefore the charge density Q v is also equal to 
the volumetric charge density of the diffuse layer [Revil 
and Florsch, 2010], 

G~ = (i-/)Gr, (ioi) 

where / is the fraction of counterions in the Stern layer and 
Qv denotes the total charge density in the pore space 
including the contribution of the Stern layer. This total 
charge density (Stern plus diffuse layers) is related to the 
CEC (in mol kg -1 ) by [Waxman and Smits, 1968], 

Qv = ftfc^JCEC, (102) 

where p s denotes the mass density of the grains. The CEC 
is another fundamental parameter describing the electro- 
chemical properties of the porous material, more precisely 
the amount of active surface sites on the mineral surface at 



a given pH value. In SI units, the CEC is expressed in 
C kg -1 , but is classically expressed in meq g -1 (with 1 
meq= 1 mmol equivalent charge, e.g., 1 x 10~ e N, where 
e= 1.6x 10~ 19 C and N is the Avogadro number, 6.022 x 
10 23 mol" 1 , 1 meq g"' = 96,320 C kg -1 ). The fraction of 
counterions / can be determined from electrical double- 
layer theory [see Revil and Florsch, 2010]. For materials 
with different types of clay minerals, the average CEC is 
determined from the respective exchange capacities of the 
constituent clay types using [Rabaute et al., 2003 ; Wood- 
ruff and Revil, 201 1], 

CEC = XkCECk + XrCEC, + x S CEC s , (103) 

where Xi represents the mass fraction of mineral i in the po- 
rous material and K, I, and S stand for kaolinite, illite, and 
smectite, respectively, in the clayey material. The CEC of 
the clay members are well tabulated (see Figure 4). If there 
is only one type of clay mineral, the CEC is given by 
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Cation exchange capacity CEC [meq g~ 1 ] 

Figure 4. Specific surface area of clay minerals S s (in m 2 
g~ ) as a function of the CEC (in meq g _1 with 1 meq 
g" '=96,320 C kg" 1 in SI units) for various clay minerals. 
The ratio between the CEC and the specific surface area 
gives the equivalent total surface charge density of the min- 
eral surface. The shaded circles correspond to generalized 
regions for kaolinite, illite, and smectite. Figure adapted 
from Revil and Leroy [2004]. The two lines corresponds to 
one to three elementary charges per unit surface area. Data 
for the clay end members are from: Patchett [1975], Lipsi- 
cas [1984], Zundel and Siffert [1985], Lockhart [1980], 
Sinitsyn et al. [2000], Avena and De Pauli [1998], Shain- 
berg et al. [1988], Su et al. [2000], and Ma and Eggleton 
[1999]. Saprolite data: Revil et al. (2013). Soil data: Chit- 
toori and Puppala [201 1]. 



CEC = ip w CECc, where <p w denotes the clay fraction (in 
weight) of the porous material and CECc denotes the CEC 
of the clay minerals present in the material. 

[43] The data set of Vinegar and Waxman [1984] data- 
base was analyzed using the linear conductivity model 
described above. The results (not shown here) are actually 
pretty close to the results of the differential effective 
medium model used by Revil [2012, Table 1]. The factor 
A(4>,m) is typically comprised between 4 and 14. In Fig- 
ure 5, the high-frequency charge density determined from 
surface conductivity is plotted as a function of the total 
charge density estimated from the measured CEC using a 
titration method. From this graph, the fraction of counter- 
ions in the Stern layer is comprised between 0.85 (85%) 
and 0.99 (99%). According to Revil [2012], the maximum 
partition coefficients are reached at high salinities with 
/(kaolinite) = 0.98, /(illite) = 0.90, and ./(smectite) = 
0.85. Revil [2012] provides a way to compute the salinity 
dependence of / using a simplified complexation model at 
the surface of the minerals. 



Ql 



O Revil et al. [2013], saprolite (rich in illite) 
• Vinegar and Waxman [1984], shaly sands 

Model Q;=(\-f)Q v 




Total charge density (C/m ) 

Qv 



CEC 



Figure 5. High-frequency volumetric charge density ver- 
sus the total charge density. The high-frequency volumetric 
charge density is determined from electrical conductivity 
measurements at various salinities (from the surface con- 
ductivity and the formation factor) while the total charge 
density is determined from the porosity and CEC. 

[44] According to Figure 5, Q v is in the range 10 5 to 10 7 
C m for permeability in the range 10" lo to 10 _1 "m z [see 
Vinegar and Waxman, 1984; Revil, 2012]. For this perme- 
ability interval and according to Figure 3, the volumetric 
charge density Q v would be in the range 1-1000 C m -3 . 
Therefore, for most porous media the assumption Q v 3> 
Q v (with the exception of shales or formations/soils 
extremely rich in clays) is justified. 



6. Comparison With Experimental Data 

[45] There is lack of experimental data to evaluate the 
accuracy and predictive power of our model with respect to 
the effect of both the water saturation and for the full range 
of frequencies covering the viscous and inertial dominated 
regimes. In the following two sections, we compare our 
model separately for (1) the frequency dependence of the 
streaming potential coupling coefficient showing how sa- 
linity and frequency affects this fundamental coupling pa- 
rameter and (2) the effect of saturation on the low- 
frequency coupling coefficient. 

6.1. Saturated Case : Effect of the Inertial Term 

[46] We first investigate the effect of salinity upon the 
streaming potential coupling coefficient. From equations 
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(27) and (38) at saturation (s w = 1) and at low frequencies 
{ujTk« 1), we obtain, 



Co 



k>Q°y 



(104) 



We first fit the values of the static coupling coefficient 
displayed in Figure 6 for the Berea sandstone using equa- 
tion (104)-(4). We obtain a s = (1.2±0.3)x 10" 3 S m" 1 and 



1.4±0.2 C m" 



These two values can be independ- 
ently confirmed using our model: (1) The value of Q° v can 
be independently obtained by equation (98), which yields 

Q v = 2.0 C m . (2) The surface conductivity can be com- 
pared to the estimate made by Moore et al. [2004] using 
electrical conductivity measurements. They found as = 
2.7 x 10~ 3 S m _1 . Therefore, there is a fair agreement 
between the present theory and the published experimental 
data. From the surface conductivity (as = 1.2x 10~ 3 S m _1 ), 
the formation factor (F= 18), and the value of the mobility 
of the counterions in the diffuse layer (/3(Na + , 25° C) = 
5.2xl0 -8 m 2 s _1 V -1 ), we can estimate the value of the 
high-frequency charge density Q v using equation (100). We 
obtain = 4x 10 5 C m" 3 . We check therefore that » 

Q v in the case of the Berea sandstone (porosity 0.23, perme- 
ability 450 mD, NaCl). This is expected because the Berea 
sandstone is a sandstone with pretty large pores (6-9 urn) 
and therefore the double layer is very thin with respect to the 
size of the pores. 
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Figure 6. Static streaming potential coupling coefficient. 
The black circles correspond to measurements by Zhu and 
Toksdz [2012] (Berea sandstone, porosity 0.23, permeability 
450 mD, NaCl). The crosses come from laboratory meas- 
urements by Moore et al. [2004] (Berea sandstone, porosity 
0.19, water). The plain line corresponds to the fit the pro- 
posed model to the data of Zhu and Toksdz [2012] only. 



[47] We used now the previous results to predict the fre- 
quency dependence of the streaming potential coupling 
coefficient of the Berea sandstone. In Figure 7, we compare 
the prediction of our model with the recent measurements 
of the dynamic streaming potential coupling coefficient 
from Zhu and Toksdz [2012] (the value reported at 1 kHz 
are actually the static values). The present model is able to 
reproduce these data very well up to 100 kHz for five dif- 
ferent salinities. 

6.2. Unsaturated and Quasi-Static Case 

[48] We can now evaluate the effect of water saturation 
upon the streaming potential coupling coefficient by substi- 
tuting inside equations (27) and (38) the volumetric charge 
density, the permeability, and the electrical conductivity by 
their expressions as a function of saturation. At high salin- 
ity, a Q ~ s n w a w lF, the low-frequency charge density scales 
as Q v /s w , and the relative permeability scaling with the 
saturation is given by equation (68). We therefore obtain 
the following expression for the quasi-static steaming 
potential coupling coefficient: The quasi-static streaming 
potential coupling coefficient is therefore given by 
Co = C r Cs, where the quasi-static streaming potential cou- 
pling coefficient at saturation is given by equation (15) and 
the relative coupling coefficient is given by, 



C,- = V , (105) 
Because the streaming potential coupling coefficient needs 
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Figure 7. Dynamic streaming potential coupling coeffi- 
cient (the value reported at 1 kHz are the static values). The 
data are from Zhu and Toksdz [2012] for the same Berea 
sandstone (porosity 0.23, permeability 450 mD, NaCl). The 
relaxation is due to the transition between the viscous-lami- 
nar flow regime at low frequency and the inertial laminar 
flow regime at high frequencies. Below 10 kHz, the stream- 
ing potential coupling coefficient can be considered inde- 
pendent on the frequency. 
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to be equal to zero at the irreducible water saturation, we 
can replace the water saturation by the effective water satu- 
ration s e . 

C r = s/, (106) 

where s/ denotes the reduced water saturation and 
r = 2/A + 2— n. The same approach at low salinity (sur- 
face conductivity dominates) yields r = 2/X+3— n. In 
Revil [2013, Appendix B], it is shown that (2 + 3A)/A 
scales as (n+2). So at high salinity, this yields r = 1. In 
Figure 8, we plot a number of recently reported experimen- 
tal data, and we confirm that indeed r = 1 can be used as a 
good first-order approximation. At low salinity, a similar 
analysis yields r=2. 

7. Application to a Remediation Problem 

[49] We show an application of the presented model to 
the detection of the NAPL(oil) water encroachment front 
during the remediation of an NAPL(oil) contaminated aqui- 
fer by flooding the aquifer with water (such remediation 
can be enhanced also with the use of surfactants, e.g., Mer- 
rier and Cohen [1990], Pope and Wade [1995], Londergan 
etal. [2001]). We will use subscript "o" to characterize the 
properties of oil. 

7.1. Simulation of Water Flooding of an NAPL- 
Contaminated Aquifer 

[50] We follow the following two steps to simulate the 
remediation of a NAPL(oil)-contaminated aquifer. 

/ "\ 

+ Revil et al. [2011] (clayrock) 

0 Vinogradov and Jackson [2012] (St Bees sandstone) 

1 Guichet et al. [2003] (Fontainebleau sand, Argon) 




0 0.2 0.4 0.6 0.8 1.0 



Water saturation, Sw (-) 

Figure 8. Quasi-static relative streaming potential cou- 
pling coefficient as a function of the water saturation. The 
plain line corresponds to the equation C r = s e with an irre- 
ducible water saturation s r =0.2. Data from Revil et al. 
[2012], Vinogradov and Jackson [2011], and Guichet et al. 
[2003]. 



[51] Step 1. We used the approach developed in Karaou- 
lis et al. [2012] to generate a 2-D heterogeneous aquifer in 
terms of porosity and permeability (Figure 9). A random 
field for the clay content was generated with the SGeMS 
library [Stanford University]. We used an isotropic semi- 
variogram to compute the clay content distribution (see 
Karaoulis et al. [2012, Appendix A]). The porosity and 
permeability were then computed according to the petro- 
physical model defined by Revil and Cathles [1999]. This 
heterogeneous aquifer is assumed to be initially saturated 
with 75% of light motor oil, resulting from an oil spill. The 
initial water saturation of the aquifer is therefore s w = 0.25, 
which corresponds to the irreducible water saturation s r . 
This aquifer is located between two wells: Wells A and B. 
Well B is located 250 m away from Well A. The reference 
position, 0(— 80,30) for the coordinate system is located at 
the upper left corner of this domain. 

[52] Step 2. Water flooding of this aquifer is simulated 
in 2.5D by injecting water in Well A (constant injection 
rate) and removing the NAPL(oil) in Well B (constant 
pressure condition). The computations are done in two- 
phase flow conditions following the same equations as in 
Karaoulis et al. [2012, Appendix A]. The properties of the 
NAPL(oil) and water are reported in Table 3. We use a 
relatively low viscosity for the NAPl(oil) as usually the 
injected water is heated to decrease the NAPL(oil) viscos- 
ity. Also the NAPL is assumed to be the nonwetting phase 
in this numerical experiment. This is realistic only just af- 
ter an oil spill for instance. Indeed, after a period of few 
years, the wettability (surface tension) of the oil can 
change because bacteria produced, for instance, biopoly- 
mers bridging the oil molecules to the surface of the grains. 
This effect is not accounted for here. After the simulations, 
we display six snapshots (T1-T6) of the oil and water satu- 
rations (Figure 10). 

7.2. Simulation of the SE Problem 

[53] For each of the snapshots shown in Figure 10, we 
simulated a SE acquisition between the two wells. The sim- 
ulation of the SE problem is also done in two phases as 
described below. Because this problem is formally different 
from the unsaturated case discussed above (two-phase flow 
problem versus unsaturated problem), we need to accom- 
modate somehow this issue as discussed below. 

[54] Step 1 . This step concerns the modeling of the prop- 
agation of the seismic waves between the two wells. We 
use material properties values given in Table 3 for the com- 
putation of the seismic properties. The bulk modulus of the 
fluid is related to the NAPL(oil) saturation by using Wood 
formula as discussed in section 3.1. above [e.g., Rubino 
and Holliger, [2012] : 

1 + ^ (1Q7) 

Kf K„ K w 

where K a and K w denote the bulk moduli of NAPL(oil) and 
water, respectively. The shear modulus is independent on the 
saturation because none of the two fluids sustain shear stresses. 
The bulk modulus of the fluid is given by equation (107), and 
the density and viscosity of the fluid is given by, 
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Figure 9. Porosity and permeability map of an NAPL (oil) contaminated aquifer between two wells 
for a water- flood simulation. 



Pf= 0 - s >v)ft> + s wP w! 



(108) 



(109) 



where p a and p, v denote the density of the NAPL(oil) and 
water, respectively, and rj 0 and r] w denote the dynamic vis- 
cosity of oil and water, respectively. The difference of fluid 
pressure between the two phases is controlled by the capil- 
lary pressure curve, which is given by the same capillary 
pressure curve used to simulate the two-phase flow problem 
(see section 7.1 and Karaoulis et al. [2012, Appendix A]). 

[55] The geometry of the model used for the computation 
of the seismic waves is shown in Figure 11. The seismic 
sources is an explosive-like source located at position So in 
Well A (Figure 11). The receivers comprise 28 pairs of 
seismic stations and electrodes (noted as Crl-Cr28), which 



are located in Well B. The separation between these 
receivers is equal to 4 m. 

[56] We first solve the poroelastodynamic wave equa- 
tions in the frequency domain, taking into account the 



Table 3. Material Properties Used in the SE Forward Model. 



Parameter 


Value 


Units 


Reference 


/'v 


2650 


kg rrT 3 


Mavkoetal. [1999] 


Pw 


1000 


kgm~ 3 


Mavkoetal. [1999] 


Po 


900 


kgm~ 3 


Karaoulis et al. [2012] 




36.5 


GPa 


Mavkoetal. [1999] 


Kfi 


18.2 


GPa 


Mavkoetal. [1999] 


G 


13.8 


GPa 


Mavkoetal. [1999] 


K w 


2.25 


GPa 


Jardanietal. [2010] 


Ka 


1.50 


GPa 


Charoenwongsa et al. [2010] 


'/„■ 


lxlO" 3 


Pas 


Jardanietal. [2010] 


■lh 


50xl0~ 3 


Pa.s 


Light motor oil 
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a. Reference - SnaphotTl b. 200 days - Snaphot T2 b. 400 days - Snaphot T3 




Distance (m) Distance (m) Distance (m) 

Figure 10. Six snapshots showing the evolution of the water saturation s w over time in a 150-m-thick 
NAPL(oil) contaminated aquifer. The initial water saturation in the aquifer is equal to the irreducible 
water saturation s r = 0.25 (which correspond to a NAPL saturation of 0.75). In this study, the NAPL (oil) 
is considered to be the nonwetting phase. 



variable saturation of the water phase. We use the (u, p) 
formulation of section 3.2 above (see Jardani et al. 
[2010] for going to the exact form of the partial differen- 
tial equations). We use the same approach as in Rubino 
and Holliger [2012]. We use the multiphysics modeling 
package Comsol Multiphysics 4.2a and the stationary 
parametric solver PARDISO to solve the resulting partial 
differential equations [Schenk and Gartner, 2004, 2006; 
Schenk et al., 2007, 2008]. The problem is solved as fol- 
low: (i) we first compute for the poroelastic and electric 
properties distribution for the given porosity, fluid perme- 
ability, and saturation distribution of the NAPl(oil) and 
water phases, then (ii) we solve for the displacement of 
the solid phase u and the pore fluid pressure p in the fre- 
quency domain. The solution in the time domain is com- 
puted by using an inverse-Fourier transform of the 
solution in the frequency wave number domain (see Jar- 
dani et al. [2010] and Araji et al. [2012] for details 
regarding the numerical procedure). 

[57] In the frequency domain, we use the frequency 
range 8-800 Hz since the appropriate seismic wave in this 
setting operates in this range and the associated electrical 
field occurs in the same frequency range. Then using this 
frequency range, we compute the inverse fast Fourier trans- 



form (FFT -1 ) to get the time series of the seismic displace- 
ments u x and u z , and the time series of the electric potential 
response ip. We use a rectangular mesh at least 10 times 
smaller than the smallest wavelength of the seismic waves. 
We have checked that this corresponds to the smallest 
mesh for which the solution of the partial differential equa- 
tion is mesh independent. The seismic source is located at 
position So(x s = 0 m, z s ~ 116 m) with a magnitude of 1.0 
x 10 4 N m. Regarding the source time function, we choose 
a Gaussian source time function that is generated with a 
delay time of 30 ms and a dominant frequency equal to 160 
Hz. For theoretical discussion on this type of source, see 
Araji et al. [2012] and Mahardika et al. [2012]. At the four 
external boundaries of the domain, we apply a 20 m thick 
convoluted perfect matched layer (CPML). The sensors 
located at Well B are 30 m away from the right-side CPML 
layer, and therefore the solution is not influenced by the 
PML boundary condition. This receiver arrangement 
mimics the acquisition that would be obtained with triaxial 
geophones and dipole antennas. CPML boundary condi- 
tions consist of a strip simulating the propagation of the 
seismic waves into free space without any reflections going 
back inside the domain of interest (see Jardani et al. [2010] 
for further details on the implementation of this approach). 
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Figure 11. Sketch of the domain used for the modeling. 
The domain isa350m x 350m square. The injection well 
(Well A) is located at position x = 0 m. The seismic source 
is located in this well at So (0 m, 116 m). The 71 receivers 
(Crl-Cr71) are located in the production Well B located at 
x = 250 m. The discretization of the domain comprises a fi- 
nite element mesh of 125 x 125 rectangular cells. Perfectly 
matched layer (PML) adsorbing boundary conditions are 
used at the borders of the domain to simulate the wave 
propagation problem with open boundaries. 

[58] Step 2. We compute the electrical problem using 
equation (27) with m = n = 2 and a surface conductivity 
a s = 0.01 S m _1 while the conductivity of the pore water is 
setup at 0.1 S m _1 . The charge density Q v is determined 



from equation (98) and the distribution of the permeability 
(see Figure 10). The source current density is determined 
with equations (106) and (107) and is therefore consistent 
with the data shown in Figure 8 and the solution for the dis- 
placement of the solid phase and the fluid pressure. Finally, 
the electrical potential distribution is obtained by solving a 
Poisson equation for the electrical potential. 

7.3. Results 

[59] The evolution of the seismic displacement and the 
electrical potential time series recorded at station Crl2 for 
each saturation profile (T1-T6) are shown in Figure 12. The 
seismic displacements are generated from the seismic source, 
which allows only for generation of P waves. In this case, 
with the porosity distribution displayed in Figure 9 and with 
the water saturation variations shown in Figure 10, the average 
P-wave velocity of profile T1-T6 is about 4800 m s _1 . There- 
fore, the P-wave arrivals in Well B is roughly the same at the 
snapshots T1-T6 (Figure 12). Figure 12 shows that SE conver- 
sions occur for each snapshot. These conversions always 
arrive earlier than the coseismic electrical field associated with 
the arrival of the P wave. They also arrive later and later as 
the water front progresses toward Well B. There is therefore a 
clear conversion mechanism at the NAPL(oil)/ water encroach- 
ment front for each of the five snapshots T2-T6. For snapshot 
Tl, since there is no saturation contrast, we do not see any 
strong SE conversions. That said, there are still some small SE 
conversions taking place at the heterogeneities in the aquifer. 

[60] Taking the saturation profile T4 into the model, we 
show that both SE conversion generated at the NAPL(oil) 
water encroachment front and coseismic (CS) electrical 
signals are shown in all receiving stations Crl-Cr28 
(Figure 13). Figure 14 shows the seismic displacement and 
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Figure 12. Evolution of the seismic displacement and the associated electric potential time series from 
receiver point Crl2 due to changes in the position of the oil-water encroachment front during snapshots 
T1-T6. The arrival of the P wave is well identified in the seismic time series. SE denotes the SE conver- 
sions occurring at any electrical and mechanical heterogeneity in the aquifer (especially at the oil-water 
encroachment front) while CS denotes the coseismic electrical field associated with the P wave. 
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Figure 13. Seismogram and electrograms from implementing saturation profile T4. (a) The seismo- 
grams reconstructed by the geophones show the P-wave propagation from the seismic source in Well A 
to recording Well B. (b) The electrograms show the coseismic electrical potential field associated with 
the P-wave (CS) and the SE conversions with a smaller amplitude and the same time arrival. 
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Figure 14. Seismogram and electrogram at receiver Crl2 (see Figure 12). Here t 0 denoted the time of 
source ignition or the delay time (30 ms). The arrival times of the SE signals and the coseismic disturb- 
ance associated with the P wave are denoted as t\ and t 2 , respectively. The strongest signal on the electro- 
gram corresponds to the coseismic disturbance associated with the P-wave propagation (see Figure 13). 
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the electrical time series for station Crl2 with information 
on the delay time t 0 , SE conversion at time t\, and the simi- 
larities of coseismic and P-wave arrival time t 2 . 

[6i] In terms of amplitudes, the type of signal measured 
here is easily recordable in the laboratory and in the field 
through stacking. Dupuis, Butler, and colleagues have 
developed methods to improve the signal-to-noise ratio in 
SE investigations, and they demonstrated that there are no 
serious issues in measuring SE conversions in field condi- 
tions [Dupuis and Butler, 2006; Dupuis et al., 2007, 2009]. 
In conclusion, we see that our numerical model implies that 
the NAPL(oil)/water encroachment front can be detected 
through SE measurements. 



m -3 ), u/(x) — iis(x) denotes the local instantaneous veloc- 
ity of the pore water with respect to the solid phase (in m 
s _1 ), x is a local position in the pore space of the material, 
and dr is an elementary volume around point M(x). Equa- 
tion (Al) is valid whatever the size of the diffuse layer with 
respect to the size of the pores. The macroscopic source 
current density Js (called the streaming current density) is 
related to the local current density j s defined with respect 
to the deformation of the solid phase by, 



J s = ^(p(x)(ii / (x) - u s (x)), 



(A3) 



(A4) 



7. Concluding Statements 

[62] The following conclusions have been reached : 

[63] (1) We have developed the first SE/electroseismic 
model in unsaturated conditions, which is valid whatever 
the thickness of the double layer with respect to the size of 
the pores. This model is based on the recently introduced 
concept of effective charge density that can be dragged by 
the flow of the pore water. This effective charge density 
can be related directly to the permeability, thereby avoiding 
the introduction of parameters that can be difficult to assess 
independently for field applications. 

[64] (2) The model has been tested with respect to exper- 
imental data as a function of the frequency at full saturation 
and as a function of the saturation at low frequencies. There 
is a need to test the model at partial saturation in the full 
range of frequencies covering the viscous-laminar and iner- 
tial laminar flow regimes. 

[6s] (3) The model can be easily modified to account for 
two-phase flow conditions for applications to remediation 
problems of NAPL/DNAPL plumes. We described a nu- 
merical model related to an application that is relevant of 
the early NAPL(oil) contamination of a clayey sand aqui- 
fer. We have shown that the oil water encroachment front 
can be remotely detected using cross-well SE conversions. 
Therefore, this opens the door to a broad range of applica- 
tions to monitor and image remotely changes in the water 
saturation in the vadose zone and the monitoring of NAPLs 
and DNAPLs plumes. 



Appendix A: Frequency Dependence of the 
Volumetric Charge Density 

[66] We first define the effective (moveable) charge den- 
sity in quasi-static flow conditions. We write the macroscopic 
charge density (charge per unit pore volume, jp C m -3 ) that 
is dragged by the flow of the pore water as Q v . This charge 
density is given by, 

Qy(u f (x) - ii s (x)) = (p(i) (*,(*) - u s (x)>, (Al) 



QdT, 



(A2) 



3s = Q v (f>(a f {x) - iis(x)), 



Js 



(A5) 



(A6) 



This is showing how the charge density Q v is not a static 
charge density but a dynamic one associated with the aver- 
age of the pore water velocity. This is why this effective 
charge density is very sensitive to the permeability of the 
porous material as shown in Figure 3. 

[67] We continue now with the thin-double layer theory 
of Pride. In this limit, our model has to be compatible with 
his theory based on the zeta potential. Taking the low-fre- 
quency limit of Pride [1994, equation (237)] and dropping 
the second-order terms in the thin double layer limit, the 
cross-coupling term L* (uj) is given by, 



L*(oS) 



sj 1 — iu)T k ' 



(A7) 



and in the present context, the low-frequency coupling 
term is given in saturated conditions by, 



(A8) 



In the thin double-layer theory based on the excess charge 
density, the frequency-dependent cross-coupling term 
L*(of) is related to the frequency-dependent permeability 
and charge density by (in saturated conditions), 



L*(u) 



k*{u)Q v {u) 



(A9) 



In addition, the frequency dependence of the permeability 
has already been established and is given by, 



1 — iuJTk 



(A10) 



From equations (A7)-(A10), saturated conditions, in the 
thin double-layer limit and in the low-frequency limit, the 
frequency dependent volumetric charge density is given by 



where the brackets denote a pore volume averaging, p(x) 
denotes the local charge density in the pore space (in C 



Q] 



lU)T k . 



(All) 
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Now in the high-frequency limit, the limit of the volumetric 
charge density Q v {uj) is given by, 

lim Q* r (u)=Qy. (A12) 

UJ»l/T t 

A simple way to connect simply and smoothly equations 
(Al 1) and (A12) is to use the following function, 

^Vr= l-A— /i 1 ■ ■ ( ai3) 

[68] Equation (A13) has been established in the thin dou- 
ble-layer limit. For the thick double-layer case, Q v — Q v 
(the same density of charges is dragged along with the pore 
water flow independently on the frequency) and 
Q v (lo) — Q v . In the thick double-layer case, there is no 
frequency dependence of the volumetric charge density as 
discussed in detail in Revil and Linde [2006]. Equation 
(A13) provides therefore the correct solution in this case 
too. In the Jjhin double-layer approximation (see Figure 2b), 
Qy >> Qv ( see examples in Jougnot et al. [2012] and 
section 5 of the main text). In the Jhick double-layer 
approximation (see Figure 2a), Q v « Q v « (1 —f)Qv (all 
the counterions of the diffuse layer are dragged by the flow 
whatever the frequency) and therefore the frequency- 
dependent excess charge density Q v is given by, 

Q v {lo,s w ) » ^ = (1 -f)Q v s- x . (A14) 
s w 

The charge density Qy can be determined from the clay 
composition, the CEC of the different clay minerals, and 
their mass fraction (see section 5.2 of the main text). 
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